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ABSTRACT 

Model surface brightness profiles based on Fokker-Planck simulations have been used 
to assess the high resolution surface brightness profile of the globular cluster NGC 
6397 of Lauzeral et al. (1992). The profile is well fitted by the model in the maximum- 
expansion phase of a gravothermal oscillation with a core radius of 0.06 pc (6"), but 
an unresolved core cannot be ruled out. A core this size is also expected in highly 
evolved clusters including significant numbers of primordial binaries. A more massive 
cluster will provide better statistics, as would star counts. The greater distance of the 
other globular clusters and the requirements of star counting will demand very high 
resolution observations to definitively resolve cores this size. 
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1 INTRODUCTION 

In globular clusters which have passed through core collapse, 
the central density profile takes the form of a power-law. At 
some point the profile must flatten out, forming a core of ap- 
proximately constant density, but whether this flattening is 
observable will depend on the underlying stellar distribution 
function and the sampling of this distribution function by 
the finite number of stars. From a theoretical perspective, 
the size of the post-collapse core of a cluster will depend 
on the nature of the heating mechanism which drives the 
re-expansion. Models which depend on binaries resulting 
from three-body encounters (Lee 1987) require much higher 
central densities before contraction is reversed than do mod- 
els using tidally-formed binaries (Statler, Ostriker & Cohn, 
1987) or with primordial binaries (Gao etal. 1991). Further, 
if gravothermal oscillations, which are observed in statistical 
numerical models such as gas sphere or Fokker-Planck codes, 
also occur in nature, then the post-collapse core would be 
much larger than is expected from a steady expansion. In 
the case of gravothermal oscillations the core radius is of 
order 1% of the half-mass radius of the cluster at maximum 
expansion (Murphy, Cohn & Hut 1990); from a few hun- 
dredths to a tenth of a parsec. A core this large may be 
observable under some conditions and its detection or non- 
detection would help in guiding further theoretical work. 
The finite number of stars in the region of interest places 
a fundamental restriction on the precision with which the 
comparisons can be made. 

Lauer et al. (1991) analyzed a (7-band HST image of the 
high-concentration cluster M15. After removing the bright 
stars, they claimed to have resolved a 2"2 core in the residual 
light profile. Grabhorn etal. (1992) fitted a Fokker-Planck 
model to the HST profile and stated that the observations 



were best matched by the model when it was in the state 
of maximum expansion. However, further consideration of 
the relationship between the models and the observations 
by Grabhorn etal. (1993) has demonstrated that for these 
M15 data it is in practice impossible to tell the difference be- 
tween the two extreme phases of gravothermal oscillations. 
This is because there are so few stars in the region of in- 
terest that sampling variations mask the difference. Monte 
Carlo representations of the light profile drawn from the dis- 
tribution functions of the two model phases overlap to the 
extent that stellar distributions drawn from the maximum- 
expansion distribution can look cuspy, and vice versa. 

Yanny et al. (1993) have done photometry on new HST 
V and / observations of the core of M15. They find that the 
star counts are consistent with both a power-law cusp and 
with a small core of radius less than l"5 (0.09 pc). They 
showed that given the large wings of the HST point spread 
function, the residual light profile is an unreliable guide to 
the cluster's structure. A nearer cluster is required to resolve 
the question of the observability of such a cluster's core. 

The recent high resolution observation by Lauzeral et al. 
(1992; hereafter LOAM) of the surface brightness profile of 
the nearby cluster NGC 6397 goes some way towards clari- 
fying this issue. NGC 6397 is in many ways the key cluster 
in this discussion. First, it has a central cusp in its surface 
brightness and surface density. This is generally taken as 
signifying that a cluster has passed through core collapse. 
Second, at the distance of NGC 6397 (2.2 kpc; Drukier 
etal. 1993), an arc second corresponds to 0.01 pc, and a 
maximum-expansion-phase core should be easily resolvable. 
Third, there is an extensive set of other observations avail- 
able to help constrain models of the cluster (Fahlman et al. 
1989, Drukier etal. 1993). 
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One of the problems which have arisen in understanding 
the cores of cusp-profile clusters is the origin of the colour 
gradients observed in them. Djorgovski et al. (1991) have ar- 
gued that the colour gradients arc due to a deficiency of red 
giants in the cluster cores. It is this hypothesis that LOAM 
have investigated with their new observations. These are 
based on data taken with the Danish 1.5 m telescope at ESO 
under conditions of 0!'95 seeing. What LOAM demonstrate 
is that the colour gradient in NGC 6397 can be eliminated 
if the light from the brighter red giants and from the blue 
stragglers is removed. The resulting surface brightness pro- 
file then loses its power-law shape in the central few arc 
seconds and flattens out. In the flnal proflle LOAM find the 
core radius to be 6". With the removal of the brightest stars 
the surface brightness is dominated by the less evolved stars 
and so should better reflect the underlying density distribu- 
tion. 

The comparison of these data with a Fokker-Planck 
model which includes gravothcrmal oscillations, will both 
investigate whether NGC 6397 is in a mgiximum expansion 
phase (as is expected) and address the question of the lim- 
its placed on the observability of post-collapse cores by the 
flnite number of stars. A distinction should be made be- 
tween a surface brightness profile, which is derived from the 
integrated light of the individual stars, and a surface den- 
sity profile, which comes from star counts. When a surface 
brightness profile is to be compared with a model, the mass- 
luminosity relation is needed to convert between the ob- 
served surface brightness profile and the model profile which 
is expressed in terms of the number of stars per unit area. 

2 MODEL COMPARISON 

As part of a project to compare Fokker-Planck models 
with observations of globular clusters (Drukier 1993) , I dis- 
covered one model in particular for NGC 6397 which gave a 
good match to mass functions at throe radii as well as the 
surface density profile of the cluster's bright stars. Due to 
the crowding of stellar images in the core of the cluster a 
relatively bright magnitude cutoff had to be used (Drukier 
et al. 1993), hence the central region of the observed surface 
density profile has too few stars to resolve anything within 
the central lO". The LOAM surface brightness profile allows 
the extension of this comparison into the very centre of the 
cluster. 

The models are based on the numerical integration of 
the isotropic, orbit-averaged, Fokker-Planck equation, as 
originally described in Cohn (1980). This is a statistical 
approach with the physical state of the cluster being given 
by a distribution function in energy space. The stellar mass 
spectrum is represented by a series of mass bins, each with 
its own distribution function. 

A tidal boundary has been applied in energy space fol- 
lowing the method of Lee & Ostriker (1987) . Given an initial 
tidal radius, the mean density of the cluster within that ra- 
dius is easily determined. At subsequent times, the tidal 
boundary in energy is given by the potential at the radius 
enclosing the same mean density. A fraction of the distri- 
bution function at energies greater than this "tidal energy" 
is removed at each time step. The fraction removed at a 
given energy is determined by the tidal stripping timescale 
and by the energy difference with respect to the tidal energy. 



The stripping rate is proportional to the cube of this energy 
difference. 

In order to drive the gravothermal oscillations, an en- 
ergy source is needed to oppose the gravitational collapse. 
In a star the energy comes from nuclear fusion. In a glob- 
ular cluster it is expected that the energy comes from the 
hardening of binary stars. As interactions with the rest of 
the stars in the cluster cause the binary to become more 
tightly bound, energy is liberated to reverse the collapse. In 
the models used here, the binaries are assumed to form via 
three-body interactions. This energy source is represented 
statistically by calculating the binary creation rate and the 
average amount of energy released by each binary. 

These models arc much the same as those in Drukier, 
Fahlman, & Richer (1992), but two additions have been 
made. First, the effects of stellar mass loss (due to stel- 
lar evolution) have been included by allowing the masses of 
the various mass bins to vary with time in much the same 
manner as in Chernoff & Weinberg (1990). This leads to a 
much more realistic handling of the mass spectrum and rem- 
nant population than when no stellar evolution is included 
(as, for example, in Drukier etal. 1992). It also restricts 
the choice of initial mass function to ones which do not re- 
sult in so much mass loss that the cluster would not survive 
to the present. This class of model is discussed further in 
Drukier (1993). At any given time, there is one mass bin 
with a changing mass. This bin, which I will refer to as 
the "evolving bin", represents the stars which currently are 
near to the turn-off or have evolved past it onto the giant 
and horizontal branches or their final states as white dwarfs 
or neutron stars. A non-evolving bin represents stars which 
are still somewhere on the main sequence or which have fin- 
ished evolving and arc now degenerate remnants. 

The other difference is that for the model described here 
the time step has been taken to be two central relaxation 
times. Such a small time step allows gravothermal oscilla- 
tions to take place. By way of comparison, a model with 
large time steps (typically several thousand central relax- 
ation times in the post-collapse phase) was also run. The 
evolution of the core radius for both models is shown in 
Fig. 1. The curve for the non-oscillating model has been 
shifted 0.16 Gyr earlier to match the time of core collapse. 
Since the details of the evolution are dependent on the initial 
conditions and the choice of time step, these models should 
not be taken as claiming a dynamical age for NGC 6397. 

The initial model contained 25 bins covering the initial 
mass range between 0.1 and 2OM0. The initial numbers of 
stars in each bin were distributed as two power-laws, one, 
for stars with M < OAMq, with mass spectral index a; = 1.5 
(where x — 1.35 for a Salpcter mass function), the other, for 
M > O.4M0, had x = 0.9 . The mass function was made 
continuous at OAMq. The stellar evolution times and final 
masses were the same as those used by Chernoff & Weinberg 
(1990) except that stellar evolution was assumed to stop at 
12 Gyr. This was done to keep the mean mass of the stars in 
the bin which was to start evolving at that time constant for 
the purposes of comparison. By doing this, the question of 
the appropriate light-to-mass ratio to use for comparison to 
the surface brightness profile need only be addressed once 
and not for each time of comparison. One limitation of 
the approach used here to handle stellar evolution is that 
while a bin is evolving it is held to contain a mixture of 
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Figure 1. Core radius as a function of time for the model 
described here. The inset shows the core radius for the whole run, 
while the main part of the figure shows the post-core-collapse 
epoch to emphasize the gravothermal oscillations. The dashed 
line corresponds to the evolution of the same model when large 
time steps, which suppress gravothermal oscillations, are used. 
The curve for the model without oscillations has been shifted 
0.16 Gyr earlier in order to match the core collapse times. The 
two vertical bars indicate the times of the two profiles shown in 
Fig. 2. 

main- sequence stars and degenerate remnants. The cutoff 
at 12 Gyr, being the assumed lifetime of a star with the 
mass of the bin boundary, removes this inconsistency from 
the comparison. On the other hand, this limit results in 
another inconsistency, one between the mass of this bin, 
which is the one contributing most of the Ught, and that 
of the stars presently observed at the turnoff and above in 
NGC 6397. The small mass difference should not have a 
large impact on the dynamics. At such late times the effects 
of stellar evolution are quite small, so the lack of the heating 
due to stellar evolution should also not be a concern. The 
initial distribution function was that appropriate to a King 
(1966) model with Wq — 6 with total mass 7.5xlO^M0 and 
limiting radius 39.1 pc. The initial tidal radius was also set 
to be 39.1 pc. 

The model was run for 10^ time steps taking about 12 

cpu-days on a SparcStation IPC. Six major re-expansions 
were observed in the 600 Myr span covered by the post- 
core-collapse phase shown in Fig. 1. In order to contrast the 
possible extremes in the size of a post-collapse core, I will 
be using the two extreme extrema indicated in Fig. 1 for the 
comparisons. 

For comparison with the observations the model den- 
sities were first projected and then integrated over annuli 
corresponding to the LOAM data points. The projected 
number densities at each radius were then converted to sur- 
face brightnesses by adopting appropriate magnitudes for 
each mass bin making a significant contribution to the clus- 
ter luminosity. There are several concerns which must be 



addressed in choosing these magnitudes. The first is the 
importance of the contribution of the low mass stars to 
the overall light of the cluster. For these low mass stars, 
logL ~ 2.5 log M + constant (Bergbusch & VandenBerg 
1992). If the local mass function has x < 1.5, then the cu- 
mulative contribution of these stars converges. If x is much 
smaller then the total light from the faint stars is unim- 
portant. In the model presented here, mass segregation, the 
enhancement of higher mass stars with respect to lower mass 
stars in the centres of clusters due to equipartition of energy, 
reduces the mass spectral index of the low mass stars from 
its initial value of a; = 1.5 to a; = 0.8 in the region we are 
interested in. Therefore we should be able to safely neglect 
the lowest mass bins. 

How many bins are actually needed was determined in 
the following manner. LOAM removed the blue stragglers 
and the stars with magnitudes brighter than V = 14.0. I es- 
timated the mean magnitude for the stars between V = 14.0 
and V = 16.4 using the luminosity function of Alcaino et al. 
(1987). This came out to be 15.5, but it is unclear whether 
the V = 16.4 magnitude cut off is appropriate for the turn- 
off bin. The reason for this is the mismatch mentioned above 
between the mass of this bin and that of the turn-off stars in 
the cluster. Further, there is a mismatch in age and metal- 
licity between the isochrone used to assign the evolutionary 
time scale of the mass bins in the cluster model and the ob- 
served age and metallicity of NGC 6397. In any event, first 
estimates for the the mean magnitudes of the next three 
mass bins were take from the isochrone which is from Van- 
denBerg and Bell (1985). 

The gravothermal oscillations in this model affect only 
the inner 0.1 pc, so I attempted to improve these magni- 
tude estimates by fitting the two model profiles to the eight 
outermost LOAM data points. These fits did not provide a 
strong constraint on the magnitudes. In the case of the least 
massive of the four bins a limit of 1^ > 18. was found. This 
suggests that there arc an insufficient number of stars in 
this bin to contribute significant light to the surface bright- 
ness profile. This mass bin has an upper mass limit of 0.58 
Mq. Given the initial mass function for the model and the 
Z = 10^* isochrone from Bergbusch & VandenBerg (1992), 
the total light contributed from stars less massive than 0.58 
Mq is only 19% of the light of aU the stars up to V = 14.0. 
Mass segregation reduces this fraction further, to less than 
12%, so only the turn-off bin and the next two were used 
in the comparisons. Based on these fits mean magnitudes 
were adopted for the the three remaining bins, with V = 16. 
being used for the bin covering the upper main sequence 
and more evolved stars. The estimated uncertainty on these 
magnitudes is 0.3. 

A surface brightness profile from a model is based on 
a continuous distribution. Assuming that the observed pro- 
file has the same underlying distribution, the observations 
represent a discrete sampling of that distribution. A proper 
comparison between the observations and the models should 
take into account the sampling uncertainties in the model 
"data points". Figure 2 presents the comparison between 
the observed surface brightness profile and those of the 
model at the two extrema. The points are the LOAM data 
and the two lines central to the hashed regions are the model 
profiles at the two extrema. The width of the hashed re- 
gions indicates the uncertainties. The closer-spaced hashing 
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Figure 2. The data of Lauzeral etal. (1992) are compared with 
the two extreme model profiles. The solid line shows the core col- 
lapse phase and the dashed line is the maximum expansion phase. 
The bounded, closely-hashed regions indicate the la uncertainties 
associated with sampling. The wider hashing, extending beyond 
the bounding lines, include the magnitude uncertainties as well. 

within the bounding lines represent Icr Poisson errors on the 

integrated model points, while the wider shadings, extend- 
ing beyond the bounding lines, include the 0.3 magnitude 
uncertainties in adopted magnitudes added in quadrature. 

To the eye, the maximum expansion profile provides a 
much better fit to the observations than does the core col- 
lapse profile. In order to quantify this comparison, I used 
Monte Carlo techniques to produce model surface brightness 
profiles consistent with Poisson statistics in the densities. 
For each mass species contributing to the surface bright- 
ness profile, and for each observing annulus, the number 
of stars of that mass "observed" in that bin was randomly 
drawn. These numbers were then used to calculate the sur- 
face brightness in that bin. This surface brightness profile 
was then compared with the parent profile and a value 
was calculated using the same uncertainties as for the ob- 
servations. The fraction of representations which had 
greater than or equal to that observed was then used to give 
the quality of the fit. 

The first column of Table 1 gives the range of LOAM 
data points (counting from the centre) used in each compar- 
ison. The full data set contains 20 points, of which the out- 
ermost eight were used in estimating the magnitudes. The 
reason for starting some comparisons at the second data 
point is coupled with the reason for the numbers listed in 
column 2. While, as was argued above, the fainter stars, 
which have been neglected in producing these profiles, do 
not contribute a significant amount of light on average, the 
Monte Carlo sampling can occasionally introduce anoma- 
lies. For the innermost point of maximum-expansion profile 
there less than three stars in each of the mass classes. Quite 
frequently, the simulated counts in this bin would be zero. 



giving an undefined surface brightness. To alleviate this 
problem, and to test the assumption that the faint stars are 
relatively unimportant, a background surface brightness for 
the centre was assumed in these cases and these are listed in 
column 2. The value of 16.89 is the observed surface bright- 
ness in the bin. The choice for a background surface bright- 
ness did not have any effect on the probabilities. Hence the 
neglect of the faint stars is justified. Alternatively, we could 
just ignore the central bin and start at the second, since 
the innermost data point only reinforces the trend already 
established by those further out. What this demonstrates 
is that the poor fit of the core-collapse profile is not exclu- 
sively due to the central datum. There is little effect except 
to make the core-collapse profiles somewhat more probable 
than the case where the central point is included; an unsur- 
prising result. 

The rightmost two columns of Table 1 repeat these com- 
parisons, but include the uncertainty in the stellar magni- 
tudes. Since the comparison is being made in terms of sur- 
face brightness, and since I have adopted the magnitudes 
based on fits to the outer part of the observed surface bright- 
ness profile, any systematic uncertainties in either the dis- 
tance to NGC 6397 or the magnitudes will balance out. On 
the other hand, the fits do allow for some tradeoff in mag- 
nitude between the three mass bins. This set of compar- 
isons allows for these uncertainties by assuming that the 
probability that a given magnitude is the correct magni- 
tude for each bin is normally distributed about my adopted 
magnitudes with a la dispersion of 0.3 magnitudes. With 
this assumption, both model profiles are more probably fit, 
but the maximum-expansion phase is still the more strongly 
favoured. The reason for the differences in the probabilities 
is that outermost points contribute much less strongly to 
the distribution when the magnitude errors are included. 
The high probabilities assigned to the maximal-expansion 
model when the magnitude uncertainty is included suggests 
that 0.3 magnitudes are a conservative estimate of these er- 
rors. The two sets of comparisons provide limiting cases on 
the magnitude uncertainties. 
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On the basis of Table 1, the core-collapsed profile cannot 
be ruled out, but the maximum expansion phase, with a core 
radius of 6" (0.06 pc), is to be preferred. 

3 DISCUSSION 

Given this probable resolution of a 6"core in NGC 6397, 
what can be said about its dynamical state? One possibility 
is that NGC 6397 had a population of primordial binaries 
which have given the cluster a large core (Gao etal. 1992). 
Alternatively, if NGC 6397 is undergoing gravothermal os- 
cillations, then it is not so surprising that we have caught it 
in the maximum expansion phase, as can be seen in Fig. 1. 

However, due to the sampling problems involved, an un- 
resolved core cannot be ruled out even for this nearby clus- 
ter. Even with the repaired HST PSF, surface brightness 
data are unlikely to help in clarifying the matter. The sur- 
face brightness profile is dominated by the brightest stars 
and there are too few of them to give a good statistical 
sample. To get around this, more stars are required. This 
could be achieved by using star counts which go further 
down the main sequence, but high resolution data will be 
needed to overcome the crowding in the core. If the dis- 
torting effects of the colour gradient can be ironed out (and 
the number of stars involved is few enough that overall star 
counts should not be affected) then a thorough set of high 
resolution star counts should unambiguously discriminate 
between the core-collapse and maximum-expansion phases 
and set firm bounds on the core radius of NGC 6397. 

There are two main points to consider in extending this 
comparison to other globular clusters. On the positive side, 
more populous clusters have more stars. This reduces the 
noise associated with sampling. On the negative side, all 
the other clusters are further away from us than is NGC 
6397, most considerably further. The greater distance and 
the higher crowding associated with higher densities increeise 
the resolution requirements for resolving the core and limits 
how far down the main sequence we can count stars. M15 
is 2.7 magnitudes brighter than NGC 6397, and all other 
things being equal, this corresponds to a improvement in 
the signal to noise of the counting by a factor of 3.5. On 
the other hand, it is 4.4 times further from us than NGC 
6397, giving an expected core radius of order l", and the 
extra gain in sampling is lost in the relatively larger area of 
the cluster we are required to average over. This core ra- 
dius is consistent with the limit set by Yanny etal. (1993). 
With a repaired HST the large wings of the current point 
spread function will disappear and the diffuse light would 
once again be available for measuring the surface brightness 
profile at high resolution. In this case discrimination on this 
basis may be possible for M15. What is required to defini- 
tively resolve these very small cores is a combination of a 
large number of stars in the core of the cluster of interest 
and high resolution. The defect of surface brightness pro- 
files is that they are hostage to the brightest stars, which 
are proportionately few in rmmber. Much better in terms of 
counting statistics would be extensive star counts probing 
the main sequence. Even higher resolution observations are 
then needed due to the high degree of crowding. And even 
with the counts in hand, there may just be too few stars 
in any given cluster. These requirements highlight the diffi- 
culty of observing post-collapse cores in all but the nearest 



clusters. 
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